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ABSTRACT 


The  Analyst' s  WorkBench  is  a  deterministic  integrated  framework  developed 
and  used  by  the  Weapons  Planning  Group  at  NAWC  China  Lake.  The  model  has 
no  stochastic  capability  which  requires  all  analysis  to  be  conducted  using 
parameters  based  on  expected  values  of  occurrence.  This  thesis  develops  a 
stochastic  enhancement  that  can  be  incorporated  in  the  Analyst's  WorkBench. 
Independent  identically  distributed  (IID)  events  can  be  generated  by  calls  to  the 
enhancement  as  a  parametric  input.  To  demonstrate  the  application  of  a  stochastic 
process  within  the  Analyst7 s  WorkBench,  a  test  scenario  of  a  ship  defense  model  is 
developed.  A  large  scale  missile  attack  is  simulated  deterministically  and 
stochastically  to  demonstrate  the  differences  of  a  random  probability  of  successful 
defense  vice  an  expected  value  of  success.  It  is  shown  that  the  stochastic  results 
provide  a  more  realistic  simulation  and  that  the  deterministic  results  overstate  the 
capability  of  a  system  subject  to  random  events  that  can  be  described  by  a 
statistical  distribution. 
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I.  INTRODUCTION 


This  chapter  is  a  review  of  the  Analyst7 s  WorkBench  and  discusses  the 
motivation  for  developing  a  stochastic  enhancement  that  could  be  incorporated 
within  the  model.  The  organization  of  the  remainder  of  the  thesis  is  presented  at 
the  end  of  the  chapter. 

A.  THE  ANALYST’S  WORKBENCH 

The  Analyst’s  WorkBench  (AWB)  is  an  integrated  framework  for  the 
interactive  application  of  models  and  tools  by  the  warfare  analyst  or  fleet  operator 
that  was  originally  developed  for  use  by  analysts  in  NAWCWPNS  Weapons 
Planning  Group  at  China  Lake  California  for  Strike,  AAW  and  ASuW  analyses.  It 
is  a  highly  flexible  and  extensible  environment  that  can  be  applied  across  a  broad 
set  of  analytical  applications.  It  allows  the  analyst  to  step  through  complex 
scenarios,  pausing  at  times,  or  events,  to  utilize  a  variety  of  analysis  tools  and 
models.  The  Analyst's  WorkBench  has  an  extensive  graphical  user  interface  which 
provides  the  user  the  ability  to  dramatically  visualize  the  analysis  as  well  as  the 
capability  to  document  analyses  directly  into  presentations  or  documents. 

The  Analyst's  WorkBench  has  been  used  primarily  at  China  Lake  in  a 
variety  of  applications.  Initially  planned  as  a  replay  tool  for  Weapons  and  Tactics 
Evaluation  Center  (WEPT AC)  games,  it  is  used  to  allow  analysts  an  easy  method 
to  interact  with  scenarios  and  play  "what-if?"  excursions.  It  has  been  used  to 
model  Integrated  Air  Defense  systems  (IADS)  to  allow  detailed  investigation  of 
strike  scenarios.  Its  most  recent  application  has  been  in  a  force  mix  study  which 
draws  heavily  on  efforts  outside  of  the  model  to  provide  parametric  inputs.  It  has 
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been  used  in  US  Navy  fleet  applications  for  quick  verification  of  strike  planning. 

The  Analyst’s  WorkBench  is  highly  modular  in  concept.  The  intent  of  the 
architecture  is  to  facilitate  analysts'  development  of  highly  specialized  models  to 
populate  a  library  for  future  analysts  to  draw  on.  The  Analyst’s  WorkBench  was 
developed  with  a  core  set  of  models  and  intended  to  be  grown  with  new  models 
as  the  uses  of  the  Analyst's  WorkBench  became  more  evident  in  the  analytical 
community.  A  model  can  be  developed  for  a  specific  or  generic  study.  The  analyst 
can  select  models  from  the  library  to  support  the  study  being  conducted.  The  core 
framework  of  the  model  contains  a  simulation  clock  that  can  be  interactively 
controlled,  an  event  queue  that  can  be  modified  on  the  fly,  a  graphic  interface  that 
can  take  DMA  data  to  provide  a  realistic  display  as  well  as  allow  spherical  or  XY 
coordinate  systems  and  a  variety  of  measurement  tools  to  give  analysts  immediate 
access  to  information  in  the  midst  of  a  simulation.  [Ref.  1] 

The  analyst s  model  is  "plugged  into"  the  Analyst's  WorkBench  framework 
by  way  of  Model  Interface  Guidelines  (MIGS).  These  MIGS  allow  the  analyst  to 
use  only  those  portions  of  the  Analyst's  WorkBench  which  are  germane  to  the 
study  resulting  in  improved  computer  utilization.  Parametric  inputs  are  provided 
through  a  scripting  function  called  AWSUM  (Analyst's  WorkBench  Script  of  User- 
defined  Movements)  and  models  can  add  commands  to  the  AWSUM  script  via  the 
MIGS  architecture.  The  advantage  of  the  MIGs  "hooks"  are  that  as  a  model  is 
developed,  the  overhead  is  significantly  reduced  because  most  of  the  common 
code  can  be  called  from  the  main  model.  Currently  the  model  contains  seven  MIGs 
"hooks"  that  provide  an  ability  for  common  events  to  be  evaluated  first  as  a 
computer  system  event,  then  as  an  Analyst's  WorkBench  event  and  finally  as  a 
model  event. 

Key  to  the  extensibility  of  the  Analyst's  WorkBench  is  the  ability  to  switch 
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modes  from  a  scripted  depiction  of  a  scenario  to  an  interactive  "what~if"  analysis. 
As  the  name  implies.  Analyst's  WorkBench  is  intended  to  be  a  tool  that  supports 
the  analyst  in  a  study  effort  and  adapt  to  the  study  needs  rather  than  scoping  a 
study  to  adapt  to  the  features  of  Analyst's  WorkBench. 

The  Analyst’s  WorkBench  is  written  as  an  Object  Oriented  Program  in 
Think™  Pascal  and  runs  in  a  Macintosh  environment.  It  can  run  on  any  Macintosh 
capable  of  32-bit  addressing  although  performance  is  obviously  enhanced  when 
run  on  68040  processors  or  on  a  PowerPC  CPU.  The  complexity  of  the  models  that 
are  plugged  into  the  framework  affect  the  amount  of  RAM  actually  required,  but 
in  general  the  Analyst's  WorkBench  will  run  reasonably  well  in  an  8  MB 
environment. 

Analyst's  WorkBench  currently  has  no  ability  to  simulate  a  stochastic 
process.  In  its  current  employment,  expected  values  are  fed  into  the  model  and 
single  runs  are  utilized  to  test  scenarios.  This  prevents  the  examination  of  the 
effects  of  randomness  of  various  components  in  systems  unless  multiple  runs  are 
conducted  and  the  input  values  are  systematically  varied.  Even  this  will  not 
adequately  capture  the  effect  of  random  events  because  the  expected  value  implies 
that  all  probabilistic  encounters  will  be  treated  equally.  Clearly  this  prevents 
examination  of  queues  in  targeting  and  radar  tracking  applications. 

The  work  of  this  thesis  is  to  add  a  stochastic  enhancement  to  the  Analyst's 
Workbench.  A  tool  set  of  common  statistic  distributions  was  written  that  can  be 
called  from  the  Analyst's  WorkBench.  To  see  the  effect  of  stochastic  modeling  on 
a  scenario  similar  to  typical  tasking  of  the  Analyst's  WorkBench,  an  Anti-Air 
Warfare/  Anti-Ship  Missile  Defense  scenario  was  built.  It  contains  a  large,  multiple 
asset  threat  and  a  defense-in-depth  system  that  is  sufficiently  stressed  to  require 
queuing  to  be  examined. 


3 


A  model  was  built  that  runs  within  the  Analyst's  WorkBench  framework.  It 
runs  multiple  runs  of  the  same  scenario  with  some  of  the  input  parameters 
randomized  to  allow  demonstration  of  the  stochastic  enhancement  that  was  added 
to  the  Analyst's  WorkBench.  Specifically,  range  and  probability  distribution 
functions  (pdf)  for  various  weapon  events  were  randomized  and  the  effect  of  this 
randomization  was  measured  against  a  deterministic  formulation. 

The  focus  of  this  thesis  is  to  describe  the  tools  built  to  allow  stochastic 
modeling  within  the  Analyst's  WorkBench  architecture  and  to  describe  the  Anti- 
Ship  Missile  Defense  Model  built  to  demonstrate  a  use  of  a  stochastic  method.  A 
typical  scenario  will  be  shown  from  both  a  stochastic  and  deterministic  standpoint 
and  compared.  From  this  effort  conclusions  will  be  drawn  as  to  the  value  of 
stochastic  modeling  within  Analyst's  WorkBench. 

B.  THESIS  ORGANIZATION 

The  remainder  of  this  thesis  is  organized  as  follows: 

1.  Chapter  II  discusses  the  process  used  to  develop  a  stochastic  plug-in,  the 
random  generator  utilized  within  the  enhancement,  and  other  statistical 
distributions  that  are  supplied  with  the  enhancement  that  can  be  used  when 
specific  distribution  of  random  variables  are  desired. 

2.  ChapterlH  discusses  the  scenario  and  the  model  that  was  developed  to 
demonstrate  the  employment  of  a  stochastic  process  within  the  simulation. 
A  baseline  scenario  and  excursions  to  the  baseline  are  presented. 

3.  Chapter  IV  is  a  presentation  of  the  results  achieved  from  running 
Analyst's  WorkBench  using  the  model  developed  for  this  thesis  including 
analysis  of  the  ability  of  the  stochastic  simulation  to  produce  statistically 
random  parameters. 
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Chapter  V  presents  the  conclusions  and  recommendations  of  the  thesis 
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II.  STOCHASTIC  PLUG  IN 


The  most  important  aspect  of  building  the  stochastic  enhancement  was  to 
build  the  statistical  package  from  which  random  variables  could  be  drawn.  Law 
and  Kelton  [Ref.  2]  provide  a  lengthy  discussion  of  the  types  of  statistical 
distributions  that  one  might  draw  from  in  order  to  translate  a  real  world  model 
into  a  theoretical  simulation.  Selected  algorithms  demonstrated  in  the  reference 
were  developed  into  Pascal  code  and  then  incorporated  as  a  unit  in  the  Analyst's 
WorkBench  model.  This  represents  the  majority  of  the  coding  effort  that  occurred 
in  this  thesis  and  is  discussed  in  detail  below. 

A.  RANDOM  NUMBER  GENERATOR 

Modelers  are  well  aware  that  random  number  generators  are  not  random  at 
all  and  in  fact  do  not  want  them  to  be  in  order  to  be  able  to  reproduce  results. 
Reference  2  presents  ample  discussion  of  the  Linear  Congruential  Generators 
(LCG)  as  a  means  to  produce  pseudo-random  numbers  for  use  in  describing 
variables.  Once  an  LCG  generates  an  output  that  appears  as  uniformly  distributed 
between  0  and  1,  additional  statistical  distributions  can  be  easily  produced. 
Clearly,  the  properties  of  the  random  number  generator  need  to  be  rapidly 
generated  variables  that  have  uniform  appearance  without  undue  correlation  and 
be  reproducible.  It  is  important  to  have  a  long  period  so  that  the  simulation  can  be 
run  long  enough  without  risk  of  repeating  numbers.  Much  discussion  is  given  to 
the  advantages  of  using  prime  modulus  multiplicative  LCGs  and  in  determining  the 
type  of  random  number  generator  to  use,  it  seemed  unnecessary  to  develop  a  new 
LCG.  In  short,  the  Marse-Roberts  code  utilized  in  Reference  2  was  sufficient  to 
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provide  very  adequate  random  number  generation. 

As  previously  stated,  the  Analyst's  WorkBench  is  a  Macintosh  computer 
based  model  and  the  code  is  written  in  Think  Pascal™  [Ref.  3].  Think  Pascal™  does 
contain  a  random  number  generator  but  in  order  to  (1)  ensure  portability  to  other 
platforms  and  (2)  ensure  adequate  control  of  the  random  number  generation,  it 
was  decided  to  utilize  the  Marse-Roberts  random  number  generator  presented  in 
Law  and  Kelton.  This  code  requires  a  32-bit  computer  architecture  to  run  properly 
and  some  older  Macintosh  platforms  are  only  16-bit  architecture.  However,  the 
Analyst's  WorkBench  is  also  based  on  the  32-bit  architecture  and  this  was 
determined  to  be  an  acceptable  solution.  The  LCG  is  described  in  Reference  2  and 
follows  the  formula: 


Zj  =  (a  ■Zi_l  + c)(modm)  2.1 

where  c= 0; 

a=  630,360,016; 

and  m=  2,147,483,647. 

The  U(0,1)  distribution  is  produced  from  Ut= Z^m.  The  Marse-Roberts  code 
provides  the  generation  of  the  Uniform  random  numbers  as  part  of  the  code  "to 
avoid  distortion  and  machine-dependent  variation  due  to  the  method  of 
computing  and  storing  floating  point  numbers."  [Ref.  2]  This  random  number 
generator  is  widely  used  and  the  empirical  chi-square,  serial  and  runs  tests  were 
well  presented  in  the  reference  so  independent  testing  to  verify  the  appearance  of 
a  Uniform  random  number  was  not  necessary.  Figure  2-1  shows  a  scatter  plot  of 
1000  random  numbers  generated  from  the  first  seed  value  in  the  code.  Visually, 
this  has  the  appearance  of  randomness  that  would  be  desired  in  a  random  number 
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The  feature  of  Marse-Roberts  code  that  made  it  particularly  appealing  was 
the  ability  to  start  100  different  streams  of  "HD  uniform  random"  variables  by 


calling  sequential  integer  seed  numbers.  Thus,  if  a  variety  of  random  variables  are 
required  as  input  parameters  within  a  simulation  replication,  determining  the  seed 
number  can  be  accomplished  by  using  an  integer  counter  rather  than  explicitly 
providing  seed  values. 

The  implementation  of  the  generating  ii(0,l)  is  accomplished  either  by 
specific  calls  within  a  specific  module  or  via  AWSUM  scripts,  when  implemented. 
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A  random  number,  X,  distributed  Uniform  (0,1)  is  called  by  rand(x)  where  x  is  the 
integer  seed  value  for  the  LCG. 

B.  STATISTICAL  DISTRIBUTIONS 

Most  simulations  employing  the  Analyst's  WorkBench  would  probably  use 
variables  from  some  well  known  statistical  distributions  where  only  the  location, 
scale  and  shape  parameters  need  to  be  defined.  In  fact,  it  is  envisaged  that  a  user 
of  the  Analyst's  WorkBench  in  a  stochastic  environment  would  either  know  in 
advance  the  type  of  statistical  distribution  the  variables  came  from  or  want  to  test 
a  hypothesis  using  several  different  distributions.  Reference  2  provides  a  fairly 
comprehensive  compendium  of  useful  probability  distributions  that  were  used  as 
the  basis  for  importing  them  into  the  stochastic  enhancement  developed  for  the 
Analyst’s  WorkBench  in  this  thesis.  The  challenge  was  to  transform  the  variables 
generated  by  the  LCG  into  random  variables  representing  these  distributions  that 
could  be  called  by  the  model.  Reference  2  suggests  that  when  generating  random 
variables  one  should  take  into  account  the  exactness  of  the  generated  variable,  the 
efficiency  of  the  generation  process  and  the  robustness  of  the  algorithm  used.  It 
discusses  different  ways  to  generate  random  variables  but  weighs  the  discussion 
heavily  in  favor  of  using  the  Inverse  Transform  Method  that  follows  the  general 
algorithm  of: 

a.  Generate  U  ~  U(0,1). 

b.  Return  X  =  F2(17). 

The  algorithms  presented  in  Chapter  Eight  of  Reference  2  were  utilized  as 
the  basis  of  the  stochastic  unit  code  presented  in  Appendix  A.  The  algorithms  as 
they  appears  in  Reference  2  are  reviewed  here  with  some  explanation  of  how  an 
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analyst  might  want  to  employ  them  with  in  the  Analyst's  WorkBench.  In  all  cases 
the  analyst  would  have  to  provide  parameters  to  bound  the  distribution;  location 
parameters,  y,  the  scale  parameters  p,  and  the  shaping  parameters  a.  These 
parameters  normally  will  be  provided  in  the  AWSUM  script  as  part  of  the  set  up 
in  lieu  of  a  specific  parameter  value  provided  to  the  Analyst's  WorkBench  as  it  is 
presently  used.  An  additional  parameter  required  is  a  seed  value.  Using  a  variety 
of  seed  values  for  the  generation  of  the  initial  uniform  random  variate  will  allow 
an  additional  if  not  arbitrary  randomization  facet. 

1.  Uniform 

The  Uniform  Distribution  is  easily  obtained  because  whenever  a  random 
number  is  called,  a  Uniform  value  distributed  between  0  and  1  is  returned.  The 
algorithm  used  to  generate  a  I7(a,b)  uniform  variate  is: 

1.  Generate  U~U(0,1) 

2.  Return  X  =  a+(b-a)U 

The  code  for  this  random  variable  is  called  by  unifomw  and  requires  the 
analyst  to  supply  a  seed  integer  for  the  random  number  generator;  and  low  and 
high  real  number  to  bracket  the  distribution. 

2.  Normal 

The  Normal  distribution  is  certainly  well  understood  and  widely  used  in  a 
variety  of  applications.  One  use  is  the  simulation  of  munitions  impact  circular 
error  probable  (CEP).  The  algorithm  to  create  the  Normal  distribution  actually 
produces  two  HD  Normal  (0,1)  random  variates  but  only  one  is  utilized  in  the  call 
in  order  to  simplify  the  scripting  call  in  the  program.  The  method  selected  to 
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generate  Normal  random  variates  is  the  polar  method  described  by  Marsaglia  and 
Bray  and  is  shown  below: [Reference  4] 

1.  Generate  Uj  and  u2  as  IID  U(0,1),  let  v,  =2^-1  for  r=l,2  and  let 

w  =  v*+v* 

2.  If  w  >  1 ,  go  back  to  step  1.  Otherwise,  let  y  =  V(-2inW)/w  ,  =  VjY,  and 

x2  =  v2y-  Then  Xj  and  X2  are  IID  N(0,1)  random  variates. 

It  is  noted  [Ref.  2]  this  method  is  slower  than  some  newer  methods,  but  the 
simplicity  of  the  method  overrides  the  speed  gains.  Of  note  with  the  Normal 
distribution  is  that  it  can  be  used  to  generate  negative  values.  This  would  not  be 
useful  in  event  queuing  calls  but  can  be  very  useful  in  placement  decisions.  If  the 
analyst  chose  to  use  a  Normal  distribution  for  event  queuing  (such  as  determining 
actual  launch  time  around  a  scheduled  Time  Of  Launch)  these  events  would  have 
to  be  scheduled  prior  to  running  the  event  queue  or  inserting  some  mechanism  to 
reject  negative  numbers.  Also,  the  IID  normal  variates  that  are  produced  do  not 
generate  a  true  normal  distribution  since  neither  the  Normal  distribution  or  its 
inverse  have  a  simple  closed  form.  However,  this  should  not  be  a  problem  in  the 
application  of  the  Analyst's  WorkBench  because  extreme  normal  variables  will 
tend  to  be  rejected  or  numerically  resolved  in  practice. 

The  analyst  must  specify  the  seed  number  for  the  LCG  and  the  location  and 
the  shape  of  the  distribution  by  supplying  a  mean  40  and  a  sigma  value  (a) .  The 
random  variable  is  then  calculated  by  X'  =  \i  +  ax.  The  code  is  called  by  Normaix. 

3.  Exponential 

The  exponential  distribution  is  very  useful  for  modeling  interarrival  times 
between  events  or  service  times  that  occur  at  a  regular  rate.  Initial  radar  detections 
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could  be  modeled  by  an  exponential  distribution  with  the  mean  describing  the 
average  time  between  detections. 

The  algorithm  for  generating  an  exponential  random  variable  is  shown 
below.  Reference  2  makes  the  point  that  the  U  used  in  step  2  is  actually  1  -U  to 
make  the  correlation  between  the  X's  and  U's  positive.  The  code  does  this  without 
action  of  the  analyst. 


1.  Generate  £/~t/(0,l) 

2.  Return  x  =  (-p)lnt/ 


The  call  to  generate  an  exponential  random  variable  is  Exponentiaix, 
specifying  the  seed  value  and  a  real  C  value. 


4.  Gamma 

The  shape  of  the  Gamma  distribution  curves  suggests  a  possible  use  to 

describe  time  to  complete  a  task.  The  scale  and  shape  parameters  ap  describe  the 
mean  time  to  complete  the  task.  There  can  never  be  negative  or  zero  values. 

Reference  2  notes  that  X'  ~  gamma(a,p)  is  merely  P(X  ~  gamma  (a,l) 
provided  P  >  0  and  that  the  gamma  (1,1)  is  the  exponential  distribution  with  mean 
1.  This  lead  to  the  algorithm  that  has  three  different  cases  dependent  on  the  value 
of  a;  0  <  a  <1,  a  =  1  and  a  >1.  In  the  case  where  0  <  a  <1,  the  algorithm  specified 
in  reference  2  as  follows: 

1.  Generate  Ui~U(0,l)  andlet  P = bU  j. If  P>l,gottostep3.  Otherwise  proceed 
to  step  2. 

2.  Let  7  =  P 1/<x,  and  generate  U2~U( 0,1).  If  U2  <  e  "y,  return  X=  Y.  Otherwise, 
go  back  to  step  I. 

3.  Let  7=  -ln[<2>  -  P)/a]  and  generate  fVt/(0,l).If  U2<Y  return  X=  Y. 
Otherwise,  go  back  to  step  1. 
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In  the  case  where  a>l.  Reference  2  describes  an  algorithm  developed  by 
Cheng  which  use  constants  that  were  incorporated  in  this  code  as  follows: 

1.  Generate  U1  and  U2  IID£/(0,1). 

2.  Let V=  alnfcy^/J  ,Y  =  aev,Z  =  U^L^andW  =  b  +  qV-Y 

3.  If  W  +  d  -  0Z  >  0 ,  return  X  =  Y.  Otherwise,  proceed  to  step  4. 

4.  If  W  >  InZ,  return  X  =  Y.  Otherwise  go  back  to  step  1. 

The  call  to  generate  a  Gamma  random  variable  is  GaamaN,  specifying  the  seed 
value  and  a  real  a  and  C  value. 

5.  Poisson 

The  Poisson  process  is  well  understood  as  a  description  of  interarrival  times 
for  most  event  queues  where  the  probability  of  the  next  event  occurring  follows  an 
exponential  distribution.  This  process  is  particularly  useful  in  modeling  the  rate 
that  events  are  placed  on  a  queue  when  the  process  leading  to  the  appearance  on 
the  queue  is  not  as  important  to  the  analysis.  The  call  to  generate  a  Poisson  random 
variable  is  poissionp,  specifying  the  seed  value  and  a  real  X  value. 

6.  Weibull 

The  Weibull  distribution  can  be  used  in  reliability  studies  to  specify  the  task 
completion  times  or  equipment  failure  times.  A  Weibull  distribution  is  useful  in 
those  instances  where  negative  values  about  the  mean  are  not  desired.  The  Weibull 
distribution  inverse-transform  algorithm: 

1.  Generate  U  IIDt/(0,l). 

2.  Return  X=  P(-lnU  )1/ot 
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The  call  to  generate  a  Weibull  random  variable  is  weibuii,  specifying  the 
seed  value  and  a  real  a  and  p  value. 

7.  w-Erlang 

Many  stochastic  processes  are  described  as  following  an  Erlang 
distribution.  The  m-Erlang  distribution  is  achieved  by  using  a  convolution 
algorithm  to  sum  a  series  of  HD  exponential  random  variables  with  the  same 
mean.  The  call  to  generate  a  m-Erlang  random  variable  is  mEriang,  specifying  the 
integer  seed  value.  Hie  integer  m,  and  a  real  P  value. 

8.  Cauchy 

The  Cauchy  distribution  is  symmetric  around  zero,  bell  shape,  but  has 
considerably  heavier  tails  than  the  normal  distribution.  The  Cauchy  distribution 
may  be  used  similarly  to  the  Normal  but  because  of  the  heavy  tails,  no  moments 
exist  in  the  Cauchy.  The  mean  of  a  collection  of  HD  Cauchy  variables  is  no  better 
an  estimator  of  the  central  tendency  of  the  distribution  than  any  of  the  individual 
observations.  The  call  to  generate  a  Cauchy  random  variable  is  cauchyx,  specifying 
the  seed  value  and  a  real  median  and  P  value. 

C.  ADDITIONAL  STOCHASTIC  PROCESSES  BUILT  FROM 
COMBINATION  OF  DISTRIBUTIONS 

Using  the  technique  of  convolution,  additional  variables  can  be  generated. 
The  easiest  method  is  to  sum  the  random  variables  from  different  distributions. 
Reference  2  provides  methodologies  for  generating  the  desired  variables  from 
combinations  of  distributions.  The  Analyst's  WorkBench  will  generate  these 
combinations  of  random  numbers  through  custom  AWSUM  script  calls. 
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III.  SCENARIO  FOR  TESTING  THE  STOCHASTIC 
ENHANCEMENT 

Once  the  stochastic  enhancement  was  added  to  the  Analyst's  WorkBench, 
an  example  scenario  was  created  to  illustrate  how  this  enhancement  might  be 
employed.  As  a  result  of  this  thesis  research,  a  stochastic  module  can  be 
implemented  in  future  revisions  of  the  Analyst’s  WorkBench  but  immediate 
implementation  of  this  capability  is  not  possible  until  supporting  changes  to  the 
stochastic  enhancements  are  provided  in  a  future  revision  to  the  Analyst's 
WorkBench.  Therefore  a  separate  testing  MGS  was  implemented  to  run  the  model 
in  a  stochastic  mode  and  to  demonstrate  that  the  stochastic  mode  would  provide 
a  different,  more  illustrative  answer  than  would  be  obtained  by  running  the 
Analyst's  WorkBench  in  a  deterministic  mode. 

A.  SHIP  DEFENSE  MODEL 

In  an  article  that  appeared  in  US  Naval  Institute  Proceedings  [Ref.  5],  a 
scenario  was  proposed  to  demonstrate  a  comparison  of  force  structure  capability 
in  the  Littoral  warfare  environment.  A  flotilla  of  small  USN  fast  attack  craft  was 
compared  with  a  task  force  of  two  USN  DDG-51  class  destroyers  to  demonstrate 
the  cost  savings  and  force  conservation  that  could  possibly  be  obtained  with  a  fast 
attack  flotilla.  In  this  scenario,  the  chief  measure  of  effectiveness  was  how  the  force 
structure  fared  in  a  concentrated  missile  attack  aimed  at  the  on  station  forces.  The 
hypothesis  was  that  a  third  world  nation  could  conceivably  launch  a  large  scale 
attack  of  20  surface  to  surface  anti-ship  missiles  against  either  a  single  DDG-51  or 
a  forward  deployed  fast  attack  force.  The  intent  of  this  scenario  was  to  show  that 
either  force  could  be  overwhelmed  and  destroyed.  This  scenario  was  the  genesis 
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for  the  model  developed  to  test  the  stochastic  enhancement  although  the  model 
was  reduce  in  scope  to  represent  a  single  ship's  defense  in  depth  against  a  fixed 
number  of  missiles. 

The  module  that  was  developed  modeled  a  theoretical  warship  with  a  tier 
defense-in-depth  system  that  would  defend  against  a  coordinated,  large  scale 
missile  attack  of  20  identical  surface  to  surface  missiles.  The  missiles  would  be 
launched  nearly  simultaneously  from  nearly  identical  distances  so  that  the 
missiles  would  arrive  at  the  ship's  defensive  perimeter  at  near  simultaneous  time. 
Each  missile  would  have  some  probability  of  acquiring  the  target  ship  and  some 
probability  of  detonating  on  the  ship.  It  was  decided  that  if  three  impacts  occurred 
on  the  ship,  the  ship  would  be  destroyed. 

The  ship  was  designed  to  have  an  outer  perimeter  search  radar  that  would 
alert  the  ship's  defensive  systems  when  an  incoming  missile  was  detected.  The 
tracking  capability  of  this  radar  was  decided  to  be  sufficient  to  prevent  it  from 
being  overwhelmed  in  this  scenario.  An  overall  battle  director  would  track  and 
prioritize  the  incoming  missiles  and  task  the  tiers  to  attack  the  incoming  missiles 
as  the  missiles  came  into  the  tier's  engagability  envelope.  The  tier's  range  of 
defense  did  not  overlap  its  adjacent  tiers.  Each  tier  had  a  finite  number  of  assets  to 
allocate  against  the  incoming  missiles  and  finite  number  of  launchers  that  could  be 
used  to  deploy  these  assets  at  a  specific  rate.  As  missiles  were  detected,  the  tiers 
would  be  tasked  to  target  and  attack  the  incoming  missile  with  an  interceptor 
missile.  As  a  launcher  fired  an  interceptor,  it  would  reload  and  fire  at  the  next 
assigned  incoming  missile.  Figure  3-1  is  a  visualization  of  this  scenario. 

This  scenario  became  an  analysis  of  what  could  be  characterized  as  an  M/G/ 
n  [Ref.  6],  Angry  Customer  Queue.  The  servers  in  each  tier  are  the  launchers  that 
would  "service"  the  incoming  "customers"  (the  missiles)  which  could  be  assumed 
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to  arrive  in  an  HD  manner.  A  manager  (OTCManager)  prioritizes  each  missile 
based  on  its  distance  to  the  ship  and  creates  a  queue  of  incoming  missiles. 


If  a  tier  could  handle  a  missile  but  could  not  assign  a  launcher  to  the  missile  at  that 
moment,  the  missile  would  remain  in  the  queue. 

If  the  missile  could  not  be  serviced  by  the  tier  before  the  missile  traversed 
the  tier's  range  of  responsibility,  the  missile  would  no  longer  be  serviced  by  that 
tier  and  would  be  assigned  to  the  next  tier.  If  the  final  tier  could  not  service  the 
missile  within  a  specific  time,  the  missile  would  leave  the  queue  and  be  scheduled 
to  impact  the  ship  (hence  the  angry  customer).  An  impact  counter  kept  track  of  the 
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number  of  impacts  and  the  model  ran  until  all  missiles  were  serviced  or  impacted 
the  ship.  If  three  missile  impacts  occurred  prior  to  all  missiles  completing  their 
flight,  the  defense  tiers  were  assumed  to  be  inactive  and  all  additional  service 
events  were  canceled.  Servicing  a  missile  modeled  reducing  the  available  assets  to 
provide  service  with  each  service  event  and  the  possibility  that  the  servicing  asset 
would  be  unsuccessful  in  servicing  (i.e.,  miss)  the  missile. 

There  was  also  the  possibility  the  missile  would  miss  the  ship  or  fail  to 
detonate  which  was  modeled  as  a  separate  probability  and  had  no  effect  on  the 
customer  queue.  This  is  important  to  the  model  because  the  measure  of 
effectiveness  for  this  model  was  the  ability  to  prevent  more  than  three  missiles 
from  detonating  on  the  ship.  A  failure  to  detonate  would  not  prevent  the  servers 
from  allocating  an  asset  against  the  missile  but  would  allow  the  model  to  continue 
to  run  its  servers. 

While  it  is  recognized  that  any  impact  against  the  ship  would  have  a 
potential  impact  on  the  performance  of  the  defenses,  this  model  did  not  attempt  to 
degrade  performance  by  impact. 

This  model  is  one  dimensional;  translating  the  missile's  positional 
information  into  a  distance  value  so  that  the  graphical  interface  became 
unnecessary.  In  fact,  the  Analyst's  WorkBench  is  not  needed  to  demonstrate  this 
model  when  run  deterministically  as  the  analysis  can  be  done  using  rudimentary 
probabilities  of  the  missile  and  ship  characteristics.  If  this  model  were  used  in  the 
Analyst's  WorkBench  environment,  positional  information  and  missile  trajectory 
would  be  more  important  to  demonstrate  the  attack  more  graphically.  The 
conceptual  nature  of  this  model  allows  a  great  deal  of  latitude  in  the  specification 
of  parameters.  This  model  could  be  more  carefully  refined  to  reflect  more  realistic 
behavior  as  desired.  None  of  the  results  in  this  model  should  be  construed  to 
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provide  definitive  answers  without  further  refinement.  However,  the  model  could 
be  an  excellent  starting  point  for  a  more  precise  study  of  a  ship's  defensive 
capabilities  using  actual  parameter  values.  In  constructing  this  model,  many  of  the 
real  world  ship  defense  characteristics  where  either  ignored  or  folded  into  the 
analysis  with  some  other  parameters  in  order  to  keep  the  scenario  simplistic 
enough  to  ease  the  validation  of  the  stochastic  module.  For  example,  electronic 
countermeasures  are  an  integral  and  significant  portion  of  a  modem  warship's 
defense  capability  but  rather  than  model  that  capability  specifically  in  this  module, 
it  was  decided  to  incorporate  that  into  the  inbound  missile's  probability  of 
acquisition. 

B.  PARAMETERS  OF  THE  MODEL 

The  parameters  that  effect  the  randomness  in  an  object' s  performance  are 
numerous  and  the  level  of  model  realism  (and  subsequent  model  complexity)  is 
improved  by  utilizing  as  many  of  these  parameters  within  the  model  as  possible. 
Sensitivity  analysis  of  those  parameters  would  provide  an  indication  of  what 
should  be  actively  modeled  and  what  should  be  incorporated  into  assumptions. 
The  intention  of  this  thesis  is  to  demonstrate  how  the  model  could  be  run 
stochastically  so  rather  than  establish  a  high  fidelity  model,  a  few  choice 
parameters  were  selected  to  establish  that  the  Analyst's  WorkBench  could  work 
with  the  stochastic  enhancement. 

1.  Missile  Parameters 

The  attacking  missiles  were  the  objects  in  this  model  that  had  the  most  stochastic 
parameters.  Each  missile  was  determined  to  be  an  HD  event  that  would  have  some 
randomness  associated  with  it.  In  the  case  of  the  missile  objects,  the  distance  from  the 
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target  and  the  time  of  the  missile  launch  were  determined  to  be  two  parameters  that 
would  impact  the  inter-arrival  times  of  the  missiles  in  the  attack  most  significantly.  Speed 
is  also  something  that  can  vary  significantly  but  as  the  missiles  were  assumed  to  be 
identical,  it  was  decided  to  hold  speed  constant.  Table  3-1  lists  the  parameters  that  were 
assigned  to  each  missile  in  the  model.  The  model  was  constructed  to  allow  changing  these 


Table  3-1  Missile  Parameters 


Number 

of 

Missiles 

Speed 

Distance  from 
target  at  launch 

Time  of  launch 

P 

A  acq 

20 

420  knots 

~N(80nm,lnm) 

~N(30secs,3secs) 

0.97 

0.83 

parameters  prior  to  run  time  but  this  feature  was  not  utilized.  Missile  parameters 
remained  constant  between  base  case  and  excursion  runs  that  will  be  discussed  later  in 
this  chapter. 

The  distance  from  a  target  can  be  a  function  of  the  navigational  accuracy  of  the 
launching  platform,  geography  of  the  region  and  man  made  impediments  such  as 
shipping  traffic.  It  is  reasonable  to  assume  that  if  a  simultaneous  attack  was  desired,  the 
position  of  each  missile  at  the  time  of  the  launch  order  would  be  equi-distant  from  the 
target.  Further,  it  is  reasonable  to  assume  that  the  impact  of  navigational  accuracy, 
geography  and  shipping  would  cause  the  position  of  launch  of  all  missiles  to  be 
distributed  normally  around  a  mean  value  selected  on  the  basis  of  the  range  of  the  missile. 
The  thesis  module  uses  a  normal  distribution  with  mean  of  80nm  and  sigma  of  1  nm  for 
this  scenario. 

Ideally,  a  simultaneous  launch  means  that  all  missiles  leave  their  launchers  at  the 
same  time.  In  practice,  there  is  some  randomness  associated  with  the  time  of  launch  that 
can  be  attributed  to  factors  such  as  electrical  circuitry,  booster  fusing  and  human  error. 
The  time  of  launch  in  this  study  was  assumed  to  be  normally  distributed  with  a  mean  h- 
hour  time  from  the  start  of  the  model  and  a  3  second  standard  deviation.  Because  a 
negative  time  number  is  possible  using  this  distribution,  the  code  was  written  to  force  all 
times  to  t  ^  0. 

The  probability  a  missile  will  acquire  the  target  is  certainly  the  function  of  myriad 
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parameters  including  missile  guidance  characteristics,  target  position  uncertainty,  flight 
profile  and  some  ship  defense  aspects  like  chaff  deployment  and  electronic 
countermeasures.  Playing  the  stochastic  aspect  of  each  of  these  factors  would  distinctly 
add  to  the  realism  but  also  adds  to  the  complexity  of  the  model.  Thus  it  was  decided  that 
a  single  Probability  of  Acquisition  value  would  be  applied  to  this  model  and  that  it  would 
considered  a  binomial  distribution  with  each  acquisition  event  run  as  a  Bernoulli  trial. 
Probability  of  Detonation  is  also  considered  as  a  Bernoulli  trial.  In  most  cases  the 
probability  of  detonation  is  based  on  the  reliability  of  the  missile.  In  the  model,  random 
uniform  numbers  are  generated  from  a  uniform  distribution  and  if  the  number  is  less  than 
the  specified  number,  the  event  is  considered  a  success. 

2.  Target  Parameters 

The  tiers  were  given  a  unique  range  of  responsibility  to  prevent  overlap  problems. 
Each  tier  had  a  finite  number  of  launchers  that  could  launch  an  interceptor  and  recycle  to 
launch  another  interceptor  provided  there  were  assets  available  in  the  tier.  Each 
interceptor  was  given  a  Probability  of  Single  Shot  Kill  ( PSSK)  that  was  played  as  a  Bernoulli 
trial.  The  interceptor's  speed  was  held  constant  within  each  tier.  Table  3-2  summarizes  the 
characteristics  of  each  tier.  The  Pssk  used  represents  an  alerted  and  defended  target. 


Table  3-2:  Parameters  of  Target  Tiers 


Tier 

Number  of 
Launchers 

Number  of 
Interceptors 

Interceptor 

Speed 

(knots) 

Max 

Range 

(nm) 

Min 

Range 

(nm 

■ 

1 

4 

20 

420 

40 

10 

0.75 

2 

2 

10 

300 

10 

3 

0.63 

3 

2 

50 

630 

3 

0.25 

0.45 

The  defense  system  utilized  a  "shoot  -  look  -  shoot"  tactic.  One  interceptor 
would  be  allocated  against  the  incoming  missile  and  no  other  action  would  be 
taken  until  the  results  of  the  Intercept  event  were  known.  It  allows  the  same  tier  to 
allocate  an  additional  asset  against  an  incoming  missile  if  an  intercept  event  can  be 
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generated  prior  to  the  attacking  missile  leaving  the  tier's  range  of  responsibility. 
This  model  assumes  instantaneous  battle  damage  assessment  (BDA)  but  could  be 
easily  modified  to  delay  or  degrade  the  knowledge  that  an  interceptor  had  missed. 

C.  MODEL  OPERATION 

This  model  was  written  to  evaluate  the  effect  of  random  variables  on  the 
proposed  scenario.  This  required  multiple  runs  of  the  scenario  and  a  statistical 
evaluation  of  the  results  of  the  data  collected.  Correlation  was  determined  to  not 
be  a  problem  in  this  particular  model  but  enough  data  points  were  desired  to 
ensure  adequate  population  of  a  histogram  for  analysis  purposes.  The  model  was 
constructed  so  that  the  analyst  could  specify  the  number  of  runs  desired  and 
change  some  of  the  scenario  parameters  to  allow  excursions  from  the  baseline  run. 
The  scenario  parameters  that  could  be  changed  for  the  attacking  missiles  (bogeys) 
were  the  position  of  launch  (mu  and  sigma),  the  time  of  launch  (mu  and  sigma), 
the  number  of  attacking  missiles  and  the  missile  speed.  The  scenario  parameters 
that  could  be  changed  for  the  target  tiers  were  the  acquisition  range  for  the  ship, 
the  maximum  and  minimum  ranges  for  each  tier,  the  number  of  interceptors  for 
each  tier  and  the  PgsK  f°r  each  tier  interceptor. 

Data  from  each  run  was  collected  and  written  to  output  files  that  were  then 
evaluated  using  data  analysis  software.  The  specific  data  that  was  of  interest  was 
the  generation  time  and  position  of  each  attacking  missile,  and  the  time  of  arrival 
at  the  ship  radar  acquisition  perimeter.  The  model  collected  data  on  the  intercept 
events  including  the  time  and  distance  of  the  events  and  the  success  or  failure  of 
those  events.  Each  run  reported  the  number  of  impacts  on  the  target  and  the 
number  of  attacking  missiles  killed.  The  final  report  of  a  scenario  run  was  the 
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average  number  of  impacts  across  all  the  runs  that  could  then  be  compared  to  a 
deterministic  calculation. 

D.  EXCURSIONS 

The  baseline  scenario  as  described  above  was  selected  based  on  a  desire  to 
show  defense  in  depth  and  the  parameters  selected  tended  to  allow  success  of  the 
model.  It  was  desired  to  also  demonstrate  the  effect  of  modification  of  parameters 
on  the  outcome  of  the  scenario.  Seven  excursion  runs  were  conducted  to  show 
degradation  of  the  ships  defensive  capability.  These  excursions  were  to 

1.  Reduce  the  Pssk  °f  each  tier's  interceptors  by  25  percent. 

2.  Reduce  the  Pssk  °f  each  tier's  interceptors  by  50  percent. 

3.  Reduce  the  engagability  range  of  each  tier  by  50  percent. 

4.  Reduce  the  engagability  range  of  each  tier  by  75  percent. 

5.  Reduce  the  engagability  range  of  each  tier  by  75  percent  and  anchor  the 
range  to  the  outer  perimeter. 

6.  Reduce  the  number  of  interceptors  in  each  tier  by  50  percent. 

7.  Reduce  the  number  of  interceptors  in  each  tier  by  30  percent. 

Of  this  list,  only  the  engagability  ranges  merit  explanation.  In  reducing  the 
enSaSability  range,  the  first  two  excursions  were  anchored  to  the  ship;  i.e.,  the 
third  tier  range  was  reduced  by  x%  and  then  starting  from  the  outer  edge  of  tier  3's 
perimeter,  the  second  tier' s  engagability  range  was  reduced  and  so  on.  However, 
these  excursions  produced  unexpected  results  so  an  additional  engagability  range 
excursion  was  run  to  start  at  the  outer  range,  reduce  its  engagability  range  by  75 
percent  and  then  allow  a  gap  in  coverage  to  exist  between  the  inner  perimeter  of 
one  tier  and  the  outer  perimeter  of  the  next  tier.  In  this  excusrion,  each  tier's  outer 
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perimeter  range  matched  the  baseline  case  but  its  coverage  range  was  reduced  by 
75  percent. 

The  baseline  case  and  the  excursions  were  each  run  100  times.  The  data 
were  collected  and  evaluated  using  a  statistical  software  package  for  the 
Macintosh.  The  average  runtime  for  each  excursion  was  approximately  30  minutes 
running  on  a  68030  processor  in  6  MBs  of  RAM.  The  data  results  will  be  presented 
in  the  following  chapter. 
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IV.  SIMULATION  RESULTS 


The  scenario  was  designed  to  demonstrate  the  difference  between  a 
deterministic  approach  and  a  stochastic  simulation  approach  to  answer  an 
analytical  question.  The  principal  measure  of  effectiveness  in  the  scenario  is  the 
number  of  attacking  missile  impacts  against  the  ship.  Three  impacts  were  deemed 
to  constitute  complete  target  destruction.  The  deterministic  results  can  be 
calculated  analytically  without  simulation  as  can  the  stochastic  results  but  with 
much  higher  calculation  complexity  and  chance  for  computational  error.  The 
stochastic  results  must  be  evaluated  to  demonstrate  the  values  generated  by  the 
model  follow  the  distribution  specified  and  to  compare  the  combined  run  results 
with  the  expected  values.  The  excursions  can  be  used  to  demonstrate  that  the 
model  is  sensitive  to  parameter  variations. 

A.  DETERMINISTIC  RESULTS 

Deterministically,  each  of  the  attacking  missiles  was  initialized  at  80NM, 
launched  at  time  0  and  allowed  to  generate  events  based  on  the  set  parameters 
specified  in  the  previous  chapter.  Thus,  applying  the  probabilities  assigned  to  the 
missiles  and  the  interceptors,  a  single  computation  of  the  expected  probability  of 
single  shot  kill  (Pssk)  required.  The  number  of  impacts  on  the  ship  is  the  sum  of 

^SSKj  ~  “  ^Tierl^ ~  ®*Tier2^^ “ ®*Tier3^  ,4.1 

Pssk  f°r  all  missiles  n  sent  against  the  ship  thus  following  a  binomial  distribution: 

20 

Impacts  =  £  PSSK  4.2 

n  =  1 
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While  these  parameters  could  easily  have  been  fed  into  the  module  and  run 
once,  the  simplistic  calculation  was  accomplished  on  a  spreadsheet. 

The  attacking  missile  Pssk  f°r  the  baseline  case  is  0.041  and  the  expected 
number  of  impacts  is  0.82.  The  variance  is  0.786.  These  are  also  the  results  for  those 
excursions  in  which  the  only  parameter  that  changes  is  the  engagement  range  for 
each  tier.  In  the  case  in  which  interceptor  inventory  is  reduced,  the  expectation 
must  be  achieved  by  stepwise  evaluation.Table  4-1  is  the  tabulation  of  results  for  a 
deterministic  evaluation  of  the  scenario. 


Table  4-1:  Deterministic  Scenario  Results 


Excursion: 

Baseline 

30% 

Interceptor 

Inventory 

Reduction 

50% 

Interceptor 

Inventory 

Reduction 

25% 

Interceptor 

Pssk 

Reduction 

50% 

Interceptor 

pssk 

Reduction 

Pssk 

0.041 

0.123 

0.160 

0.123 

0.274 

Impacts 

0.82 

2.46 

3.20 

2.46 

5.48 

These  results  indicate  that  a  20-missile-attack  as  outlined  above  would 
have  an  expected  success  as  shown  and  that  additional  missiles  added  to  the  attack 
would  have  additional  success  in  a  linear  manner.  These  numbers  can  also  be  used 
to  determine  the  probability  that  three  or  more  missiles  would  impact  the  ship  by 
computing  the  binomial  distribution  using  the  expected  PSSK  and  n  =  20  as  shown 
in  Figure  4-1  below.  This  indicates  that,  using  the  stated  parameters,  the 
probability  P(Number  of  Impacts^  3)  =  0.047.  This  probability  could  be  used  as  the 
determination  of  the  effectiveness  of  the  ship's  defensive  system. 
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Deterministic  pdf 


Figure  4-1:  Expected  pdf  of  Baseline  case 


B.  STOCHASTIC  SIMULATION  RESULTS 

This  section  discusses  the  verification  used  to  show  that  the  model  behaved 
stochastically.  It  was  desired  that  the  missile  generation  time  and  distance  should 
appear  to  be  randomly  distributed  about  a  mean  generation  time  and  distance. 
Data  was  also  collected  to  determine  the  distribution  of  the  interarrival  times  of  the 
missiles  at  the  outer  edge  of  the  defense  perimeter. 

Running  one  hundred  replications  of  the  model  produced  2000  missile 
events,  which  provided  sufficient  data  to  ensure  statistically  measurable  results. 
The  simulation  is  a  terminating  simulation  in  that  the  there  were  a  finite  number 
of  missile  events  and  a  finite  number  of  interceptors  in  each  tier.  The  number  of 
replications  were  selected  to  give  sufficient  random  events  to  give  each  attacking 
missile  ED  characteristics.  Each  excursion  was  run  with  the  same  attacking  missile 
parameters  so  validation  of  the  performance  of  the  baseline  process  was  sufficient 
for  the  entire  study.  Statistical  validation  included  standard  statistical  evaluation 
of  the  computed  parameters.  The  data  was  collected  as  part  of  the  model  output 
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and  then  evaluated  using  StatView  4.5™  [Ref.  7],  a  Macintosh  based  statistical 
analysis  package. 

The  simulation  was  designed  to  instantiate  the  attacking  missile  with  a 
random  launch  time  generated  from  a  normal  distribution  N(30, 3).  The  histogram 

in  Figure  4-2  demonstrates  the  launch  time  for  the  baseline  model. 


Histogram 

of  Bogey  Generation  Time  (in  seconds) 


Time 


Figure  4-2:  Bogey  Generation  Time  Histogram  (2000  Replications) 

The  attacking  missile's  distance  from  the  target  ship  was  distributed 
~N(80nm,  lnm).  The  histogram  in  Figure  4-3  demonstrates  the  results  of  the 
replications  and  shows  that  the  parameters  to  be  distributed  as  desired. 


Histogram 

of  Bogey  Generation  Distance  (in  NM) 


Distance 


Figure  4-3:  Bogey  Generation  Distance  Histogram  (2000  Replications) 
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The  interarrival  times  of  the  attacking  missiles  where  not  modeled 
explicitly  but  the  data  was  collected  and  is  shown  in  Figure  4-4  below. 
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Figure  4-4:  Outer  Tier  Bogey  Interarrival  Times  (2000  Replications) 

The  missile  interarrival  times  appear  to  follow  an  exponential  distribution. 
This  was  verified  by  using  a  Q-Q  plot  to  compare  against  a  random  series 


generated  from  an  exponential  distribution  with  the  same  mean  as  the  calculated 
mean  of  the  interarrival  times.  As  shown  in  Figure  4-5,  the  Q-Q  plot  is  fairly  linear 
at  the  lower  ends  indicating  that  the  interarrival  time  values  roughly  follow  the 
exponential  distribution  with  a  mean  of  2.276.  This  distribution  could  be  used  as 
the  parameter  for  determining  when  missiles  arrived  at  the  outer  tier  of  the  ship's 
defensive  system.  Using  the  interarrival  time  as  the  stochastic  parameter  would 
eliminate  the  need  to  explicatly  model  the  starting  point  of  the  missile.  However, 
most  current  analysts  are  interested  in  the  complete  flight  path  of  incoming 
missiles  so  that  the  actual  time  of  detection  can  be  measured.  This  necessitates  the 
modeling  of  the  point  of  origin  and  time  of  launch  for  an  attacking  missile  event. 
The  interarrival  times  determine  queue  loading  for  the  defensive  system  and  as 
such  are  a  useful  measure  for  analysis. 
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Outer  Tier  InterArrival  Time 


Figure  4-5  Q-Q  Plot  Comparison  of  Interarrival  Time  (2000 

Replications)  to  Exponential  Distribution  (/i  =  1/2.276) 

In  general,  the  model  behaved  as  expected  in  generating  random  variables 
from  the  appropriate  distributions. 

C.  MODEL  RESULTS 

This  section  compares  deterministic  and  stochastic  results.  First  a 
comparison  of  the  base  cases  is  presented  followed  by  an  analysis  of  excursion 
cases. 


1.  The  Base  Case 

Results  were  from  the  deterministic  case  as  shown  below.  The  results  are 

Table  4-2:  Comparison  of  Deterministic  Results  to  Stochastic  Simulation 


Mean  Number  of  Impacts 

Variance  of  Impacts 

Deterministic  Results 

0.82 

0.786 

Simulation  Results 

0.51 

0.697 
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similar  but  the  mean  number  of  impacts  is  slightly  different  and  the  variance  of  the 
simulation  appears  less  than  in  the  deterministic  case.  The  impact  histogram  for 
the  base  case  (Figure  4-6)  shows  the  relative  frequency  of  impacts  which  can  be 
interpreted  as  the  probability  distribution  for  the  base  case  under  the  stochastic 
conditions  at  initialization.These  results  indicate  that  there  is  even  less  probability 
that  any  given  20-missile-attack  will  produce  3  or  more  impacts.  In  most  cases,  the 
defensive  system  will  be  sufficient  to  prevent  any  attacks  from  hitting  the  ship  (i.e., 
P(number  of  impacts  <  1)  =  0.65)  and  most  impacts  will  be  less  than  that  required 
to  kill  the  ship. 


Experimental  Results  of  100  Model  Replications 
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Figure  4-6  pdf  Histogram  of  Stochastic  Simulation  Run  (Base  Case)  j 

The  simulation  also  provides  a  better  understanding  of  the  behavior  of  each 
tier's  interceptors  that  can  be  readily  compared  to  excursion  cases.  The  base  case 
interceptor  success  rate  is  shown  below  and  demonstrates  that  the  results  are  very 
close  to  the  parameters  specified  at  initialization.  The  overall  success  rate  for  the 
interceptors  was  70.6  percent.  This  can  be  interpreted  as  the  ship's  probability  of 
shooting  down  an  attacking  missile. 
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In  the  base  case  the  mean  number  of  Interceptors  used  was  27.46  with  a 
variance  of  7.564.  Figure  4-8  shows  the  distribution  of  interceptor  usage  by 
replication.  This  Figure  demonstrates  a  fairly  uniform  band  between  22  and  34 
missiles. 
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2.  Excursion  Analysis 

The  excursions  can  be  grouped  into  three  main  categories:  reduction  of 
interceptor  probability  of  single  shot  kill,  reduction  of  interceptor  engagement 
range  and  reduction  of  interceptor  inventory.  The  excursions  in  which  the  Pssk  °f 
the  interceptors  was  reduced  showed  the  most  significant  increase  in  ship  impacts. 
Interceptor  inventory  reduction  also  results  in  significant  changes  in  ship  impacts. 
Only  the  case  in  which  engagement  range  was  modified  remained  insensitive  to 
changes.  This  is  most  likely  an  artifact  of  modeling  the  recycle  rate  of  the  launchers 
too  rapidly.  All  attacking  missiles  had  an  opportunity  to  be  engaged  by  at  least  one 
interceptor  in  each  tier  before  the  missile  had  traversed  very  far  through  the  tier. 
Missile  engagement  was  modeled  as  a  uniform  probability  across  the  tier 
enSa8ement  range.  Had  the  detection  probability  been  modeled  as  a  function  of 
decreasing  range,  the  model  may  have  demonstrated  more  sensitivity  in  the 
reduced  range  excursions. 

a.  Reduction  of  Interceptor  Probability  of  Single  Shot  Kill 

Reducing  Pssk  °f  each  interceptor  results  an  increased  number  of 
impacts  in  both  the  deterministic  and  the  stochastic  excursions.  Table  4-3 


Table  4-3:  Reduction  of  Interceptor  Pssk  Excursion 


Excursion 

Deterministic 
Expected  Number 
of  Impacts 

Stochastic 
Simulation  Mean 
Number  of  Impacts 

Baseline  Case 

0.81 

0.51 

1/4  Reduction  of  Interceptor  PgsK 

2.46 

3.25 

1/2  Reduction  of  Interceptor  Pssk 

5.48 

7.43 
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demonstrates  that  while  the  Baseline  case  has  lower  number  of  impacts  in  the 
stochastic  case  versus  the  deterministic  case,  the  deterministic  model  is  less 
sensitive  to  variance  of  the  interceptor  PSSK  parameter  than  the  stochastic  model. 
More  PSSK  reduction  excursions  would  probably  demonstrate  that  the  number  of 
impacts  on  the  ship  would  increase  more  rapidly  in  the  stochastic  simulation. 

b.  Reduction  of  Interceptor  Engagement  Range 

In  the  deterministic  case,  altering  the  distance  over  which  an 
interceptor  is  able  to  engage  an  incoming  missile  should  have  no  effect  in  the 
expected  number  of  impacts  on  the  ship.  No  parameters  that  effect  the  probability 
of  a  missile  impact  are  changed  in  this  excursion.  Thus  in  each  case,  one  would 
expect  the  results  to  be  the  same.  The  mean  and  variance  of  each  excursion  are 
shown  in  Table  4-4.  Note  that  while  all  the  stochastic  simulation  values  are 


Table  4-4:  Range  Excursion  Statistics 


Excursion 

Mean 

Variance 

Deterministic  Case 

0.82 

0.786 

Baseline  Case 

0.51 

0.697 

Half  Engagement  Range 

0.52 

0.676 

One  Quarter  Engagement  Range 

0.60 

0.788 

One  Quarter  Engagement  Range 

0.53 

0.898 

generally  the  same,  they  are  all  less  than  the  deterministic  calculation.  This 
indicates  that  the  ship's  defensive  systems  are  generally  more  successful  in 
defeating  incoming  missiles  than  would  be  indicated  by  a  deterministic 
calculation.  The  stochastic  evaluation  of  this  excursion  is  shown  below  in  Figure  4- 
9. 
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c.  Reduction  of  Interceptor  Inventory 

Intuitively,  interceptor  inventory  is  significant  in  the  overall  success 
of  the  ship  defense  system.  If  the  ship  can  not  sustain  a  counter  attack  against  the 
incoming  missiles,  more  impacts  will  occur.  Table  4-5  shows  the  comparison  of  the 
deterministic  calculations  and  the  stochastic  simulation.  Again  the  rate  of  change 
in  the  stochastic  simulation  is  more  rapid  than  in  the  deterministic  simulation, 
indicating  a  higher  sensitivity  to  changes  in  the  interceptor  inventory.  Note  that 
the  expected  number  of  impacts  in  the  50  percent  inventory  reduction  case  is  3.2. 
In  the  simulation  results,  the  minimum  number  of  impacts  were  three,  which 
would  imply  under  the  specified  scenario  that  there  was  a  very  low  probability  the 
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ship  would  escape  destruction.  Figure  4-10  shows  the  range  of  numbers  of  impacts 
in  the  excursion  simulations. 


Table  4-5:  Reduction  of  Interceptor  Inventory  Excursion 


J  Excursion 

Deterministic 
Expected  Number 
of  Impacts 

Stochastic 
Simulation  Mean 
Number  of  Impacts 

Baseline  Case 

0.81 

0.51 

30%  Reduction  of 
Interceptor  Inventory 

2.46 

2.18 

50%  Reduction  of 
Interceptor  Inventory 

3.20 

5.70 

Figure  4-10:  Box  Plot  Comparison  of  Simulation  Results 
for  Inventory  Reduction  Excursion 
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V.  CONCLUSIONS 


The  history  of  the  Analyst7 s  WorkBench  and  its  limitations  as  a 
deterministic  model  have  been  presented  in  previous  chapters.  The  development 
of  a  plug— in  that  would  allow  the  Analyst's  WorkBench  to  produce  stochastic 
results  and  a  test  scenario  has  been  shown.  The  analysis  of  the  results  has  shown 
that  different  results  occur  within  the  scenario  depending  on  whether  the  model  is 
run  stochastically  or  deterministically.  This  chapter  provides  observations  and 
recomendations  for  implementation  and  future  research  efforts. 

A.  OBSERVATIONS 

The  Analyst7 s  WorkBench  currently  provides  a  good  framework  for 
conducting  analytical  studies  in  which  the  output  values  are  deterministic.  Explicit 
modeling  of  a  handful  of  entities  can  be  accomplished  easily  and  their  input 
parameters  statistically  varied  in  order  to  perform  sensitivity  analysis.  Large  scale 
modeling  using  changes  to  input  parameters  quickly  becomes  unwieldy  as  the 
number  of  parameters  increase.  Modeling  this  way  is  useful  for  visualization  of  a 
few  parameters  but  is  not  effective  for  evaluating  a  large  scale  operation.  The 
addition  of  the  stochastic  enhancement  allows  the  Analyst7 s  WorkBench  to  be 
modified  so  that  a  large  number  of  random  processes  can  be  evaluated  easily  and 
visualized  effectively. 

Using  the  ship  defense  scenario  and  its  excursions  to  generate  output  using 
the  Analyst's  WorkBench  in  a  deterministic  mode  can  be  compared  to  its 
performance  in  the  stochastic  mode.  In  the  deterministic  case,  for  example,  the 
only  way  to  examine  engagement  activity  would  be  to  explicitly  script  the 
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engagements.  Figure  5-1  shows  the  collection  of  all  the  engagement  distances  in 
30  T - - - r 


Observations 


Figure  5-1:  BaseLine  Excursion  -  Interceptor  Achieved  Engagement  Distance 

the  baseline  stochastic  excursion.  Deterministically,  all  engagements  would  occur 
at  the  same  engagement  distance  with  the  only  variation  occurring  as  a  function  of 
the  reload  capability  of  the  tier.  Every  engagement  should  result  in  some  damage 
to  the  missile  and  every  missile  should  achieve  some  damage  against  the  ship.  This 
would  be  the  same  regardless  of  the  excursion.  Clearly  in  the  stochastic  example, 
the  actual  engagement  range  is  distributed  across  the  allowable  engagement  range 
and  can  be  attributed  to  the  interarrival  times  of  the  missiles,  the  reload  times  of 
the  interceptors,  and  the  probability  of  missile  engagement.  The  damage  to  each 
missile  and  to  the  ship  are  discrete  events  that  are  independent  of  previous  missile 
events.  Intuitively,  this  is  a  more  realistic  simulation. 

Table  5-1  demonstrates  the  differences  between  a  purely  deterministic 
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model  and  a  stochastic  model  with  respect  to  the  ship  defense  scenario.  In  all  cases 
the  stochastic  model  can  demonstrate  each  of  the  aspects  of  distribution  in  a 
manner  the  deterministic  model  cannot.  In  some  cases  such  as  the  mean  number 
of  impacts  against  a  ship,  the  values  of  the  deterministic  calculation  will  be 
significantly  different  from  results  obtained  stochastically  and  in  others  the  model 
will  be  insufficient  to  demonstrate  a  feature  that  has  analytical  significance. 


Ability  of 

Deterministic  Model 

Mean 

Distribution  of 
Engagement 

Baseline 

✓ 

X 

Reduced  Interceptor 

✓ 

X 

Engagement  Range 

Reduced  Interceptor 

X 

X 

Inventory 

Reduced  Interceptor 

X 

X 

PSSK 

Table  5-1:  Ability  of  Deterministic  Model  vs.  Stochastic 


This  is  not  meant  to  suggest  that  the  Analyst's  Workbench  as  it  is  employed 
today  is  improperly  used.  Indeed,  as  a  game  replay  mechanism,  it  is  an  easily 
manipulated  tool.  The  ability  to  visualize  the  design  of  modeling  methodologies 
makes  a  deterministic  case  very  appealing.  The  "what  if?"  capability  currently 
provided  is  essential  for  examining  a  question  raised  as  the  game  is  replayed. 
However,  the  ability  to  adequately  examine  a  range  of  possibilities  under  a 
random  environment  is  severely  limited.  A  stochastic  analysis  better  answers 
more  specific  questions  such  as,  "How  does  changing  my  defense  tier's 
engagement  range  affect  my  defensive  capability?" 
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B.  RECOMMENDATIONS 


The  stochastic  enhancement  developed  in  this  thesis  is  easily  incorporated 
into  future  version  of  the  Analyst's  WorkBench.  MIGS  hooks  allowing  the  switch 
to  multiple  runs  of  the  same  experiment  should  be  provided  in  the  basic 
framework.  AWSUM  script  language  should  be  modified  to  allow  the  analyst  to 
specify  the  statistical  distribution  from  which  a  parameter  can  be  drawn. 

The  developers  of  Analyst  WorkBench  may  desire  a  visualization  of  the 
output  from  series  of  simulation  runs.  This  would  require  the  development  of 
more  specific  data  collection  tools  that  would  run  at  the  framework  level  of  the 
Analyst's  WorkBench.  Data  would  then  be  collected  from  each  run  of  the 
excursion.  The  statistical  means  of  the  excursion  would  then  be  used  as  a 
parameter  for  the  visualization  run.  The  final  step  of  an  excursion  would  be  the 
replay  of  the  visualization  run.  The  ability  to  look  at  the  parameter  value  as  the 
visualization  is  run  needs  to  be  modified  to  show  variance  and  range  of  the 
parameter  as  a  minimum.  Additional  software  design  is  required  to  explore  this 
possibility. 

The  scenario  developed  for  this  thesis  is  used  for  illustrative  purposes  only 
but  may  be  useful  in  exploring  the  ability  of  a  ship  to  defend  itself.  Further 
validation  of  die  ship  defense-in-depth  process  would  result  in  a  more  realistic 
model  that  would  be  useful  in  studying  defense-in-depth  as  a  stochastic  process. 
All  ship  defense  capacity  is  a  trade-off  in  ship  offensive  capability. 

A  recent  thesis  at  the  Naval  Postgraduate  School  [Ref.  8]  suggested  that  the 
employment  of  chaff  and  electronic  countermeasures  were  essential  in  increasing 
the  likelihood  of  successful  ship  defense.  This  observation  might  provide  a 
rationale  for  reducing  ship  defensive  weapon  inventory  and  allowing  more 
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offensive  weapon  inventory.  The  model  developed  for  this  thesis  combined  with 
the  scenario  of  Reference  5  would  provide  useful  trade-off  analysis  to  determine  a 
ship  weapon  inventory  and  operational  employment  to  ensure  maximum 
survivability  while  allowing  mission  capability. 
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APPENDIX.  STOCHASTICS  UNIT  CODE 


This  chapter  contains  the  Pascal  code  used  to  create  the  stochastic  enhancement  to 
the  Analyst's  WorkBench. 
unit  Stochastics; 


{The  beginning  work  for  a  stochastics  unit  to  be  incorporated  in  the  Analyst  Workbench 
Framework} 

{Created:  8  December  1993  by  LCDR  Andrew  W.  Melton,  USN} 

{  Naval  Postgraduate  School} 

{  Monterey,  CA} 

{Revision  History:  3  Aug  1994} 

{Incorporated  Inverse  Transforms  from  Law  &  Kelton  "Simulation  Modeling  &  Analysis"  Chp8, 

{29  July  1995  Added  location,  shape  and  scale  parameters  to  several  random  variables} 
{17  July  changed  to  Unit  to  be  added  to  AWB} 

interface 


uses 

Sane; 

var 

index,  nl:  Integer; 

Zset,  Tset:  Longint; 

Zmg:  array  [  1. .  100]  of  longint; 


function  fuzzyReal  (r:  REAL):  REAL;  {An  easy  test  for  model  development} 
procedure  Randdf ; 

function  UniformN  (x:  integer;  a,  b:  real):  real; 
function  Exponential  (x:  integer;  lambda:  real):  real; 
function  gamma X  (x:  integer):  extended; 

function  GamtnaN  (x:  integer;  alpha.  Beta:  real):  real; 
function  Weibull  (x:  integer;  alpha.  Beta:  real):  real; 
function  PoissonP  (x:  integer;  lambda:  real):  real; 
function  PoissonNextEvent  (lambda:  real;  x:  integer):  real; 
function  mErlang  (x,  m:  integer;  Beta:  real):  real; 

function  NormalX  (x:  integer;  mu,  sigma:  real):  real;  {Not  a  true  Normal  but  close 

enough } 

function  CauchyX  (x:  integer;  median,  beta:  real):  real; 
function  Rand  (Stream:  INTEGER):  real;  {Generate  the  next  random  number} 
procedure  Randst  (Zset:  Longint;  Stream:  INTEGER); 
function  Randgt  (Stream:  Integer):  Longint; 


implementation 


type 

ZStatsPtr  =  *ZStats; 
ZStats  =  record 

Z value:  longint; 
Zindex:  Integer; 
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next;  ZStatsPtr; 
end; 

{This  is  a  random  number  generator  taken  from  Law  &  Kelton’s  book  Simulation  Modeling 
and  Analysis, 1991} 

{It  is  a  PMMLCG  supporting  100  streams  of  uniform  random  numbers } 

procedure  Randdf;  {Ihis  may  be  something  sore  appropriately  stored  in  one  of  the  globals 
units} 

var 

i:  INTEGER; 

begin {set  the  seeds  for  all  100  streams} 


Zmg[l]  :=  1973272912; 
Zmg[4]  ;=  1280689831; 
Zmg[7]:=  913566091; 
Zmg[10]  :=  604901985; 
Zmg[13]  :=  824064364; 
Zmg[16]  ;=  75253171; 
Zmg[19]:=  233217322; 
Zmg[22]  :=  403498145; 
Zmg[25];=  762430696; 
Zmg[28]  :=  76271663; 
Zmg[31]:=  336157058; 
Zmg[34]:=  595778810; 
Zmg[37]:=  6891 1991; 
Zmg[40]  :=  622401386; 
Zmg[43]  :=  1774806513; 
Zmg[46]:=  781301 10; 
Zmg[49]  :=  1351423507; 
2mg[52]  :=  922510944; 
Zmg[55]  :=  243649545; 
Zmg[58]  :=  403188473; 
Zmg[61]:=  498067494; 
Zmg[64]  :=  597104727; 
Zmg[67]  :=  536444882; 
Zmg[70]  :=  67784357; 
Zmg[73]  :=  119025595; 


Zmg[2]  :=  281629770; 
Zmg[5]:=  2096730329; 
Zmg[8]  :=  246780520; 
Zmg[ll]  :=  1511192140; 
Zmg[14]  :=  150493284; 
Zmg[17]  :=  1964472944; 
Zmg[20]  :=  1911216000; 
Zmg[23]  :=  993232223; 
Zmg[26]  :=  1922803170; 
Zmg[29]  :=  413682397; 
Zmg[32]  :=  1432650381; 
Zmg[35] 877722890; 
Zmg[38]:=  2088367019; 
Zmg[41]:=  2122378830; 
Zmg[44]:=  2132545692; 
Zmg[47]  :=  852776735; 
Zmg[50]  :=  1645973084; 
Zmg[53]  :=  2045512870; 
Zmg[56]  :=  1004818771; 
Zmg[59]  :=  372279877; 
Zmg[62]  ;=  2087759558; 
Zmg[65]  :=  1530940798; 
Zmg[68]  :=  1663153658; 
Zmg[71]  :=  1432404475; 
Zmg[74]:=  880802310; 


Zmg[3]  :=  20006270; 
Zmg[6]  :=  1933576050; 
Zmg[9]  :=  1363774876; 
Zmg[12]  :=  1259851944; 
Zmg[15]  :=  242708531; 
Zmg[18]  :=  1202299975; 
Zmg{21]:=  726370533; 
Zmg[24]  :=  1103205531; 
Zmg[27]  :=  1385516923; 
Zmg[30]:=  726466604; 
Zmg[33]  :=  1120463904; 
Zmg[36]  :=  1046574445; 
Zmg[39]  :=  748545416; 
Zmg[42]:=  640690903; 
Zmg[45]  :=  2079249579; 
Zmg[48]:=  1187867272; 
Zmg[51]  :=  1997049139; 
Zmg[54]  :=  898585771; 
Zmg[57]  :=  773686062; 
Zmg[60]  :=  1901633463; 
Zmg[63]  :=  493157915; 
Zmg[66]  :=  1814496276; 
Zmg[69]  ;=  855503735; 
Zmg[72]:=  619691088; 
Zmg[75]  :=  176192644; 
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Zmg[76]:=  11 16780070; 
Zmg[79]  :=  1142483975; 
Zmg[82]:=  786262391; 
Zmg[85]  :=  1923011392; 
Zmg[88];=  1147297105; 
Zmg[91]  ;=  190641742; 
Zmg[94]:=  620389253; 
Zmg[97]:=  364849192; 
Zmg[100]:=  547070247; 


Zmg[77]:=  277854671; 
Zmg[80]  ;=  2026948561; 
Zmg[83]  :=  1792203830; 
Zmg[86]  :=  1433700034; 
Zmg[89]:=  539712780; 
Zmg[92]  :=  1645390429; 
Zmg[95]  :=  1502074852; 
Zrag[98]  ;=  2049576050; 


Zmg[78]  :=  1366580350 
Zmg[81]  :=  1053920743 
Zmg[84]  :=  1494667770 
Zmg[87]  :=  1244184613 
Zmg[90]  :=  1545929719 
Zmg[93]  ;=  264907697; 
Zmg[96]:=  92771 1160; 
Zmg[99]:=  638580085; 


end; {Randdf } 

function  Rand  (Stream*  INTEGER)*  real;  {Generate  the  next  random  number} 

const 

B2E15  =  32768; 

B2E16  =  65536; 

MODULUS  =  2147483647; 

MULTI 1  =  24112; 

MDLTI2  =  26143; 

var 

Hil5,  Hi31,  Lowl5,  LowPrd,  Ovflow,  Zis  LONGINT; 

begin  {Rand} 

Zi  s=  Zrng[ Stream]; 

Hil5  :=  Zi  div  B2E16; 

LowPrd  *=  (Zi  -  Hil5  *  B2E16)  *  MULTI1; 

Lowl5  *=  LowPrd  div  B2E16; 

Hi31  *=  Hil5  *  MULTI 1  +  Lowl5; 

Ovflow  :=  Hi31  div  B2E15; 

Zi  *=  (((LowPrd  -  Lowl5  *  B2E16)  -  MODULUS)  +  Hi31  -  Ovflow  *  B2E15)  *  B2E16  + 

Ovflow; 

if  Zi  <  0  then 
begin 

Zi  s=  Zi  +  MODULUS; 

end; 

Hil5  *=  Zi  div  B2E16; 

LowPrd  *=  (Zi  -  Hil5  *  B2E16)  *  MULTI2; 

Lowl5  :=  LowPrd  div  B2E16; 

Hi31  *=  Hi 15  *  MULTI2  +  Lowl5; 

Ovflow  :=  Hi31  div  B2E15; 

Zi  *=  (((LowPrd  -  Lowl5  *  B2E16)  -  MODULUS)  +  Hi31  -  Ovflow  *  B2E15)  *  B2E16  + 

Ovflow; 

if  Zi  <  0  then 
begin 

Zi  *=  Zi  +  MODULUS; 

end; 

Zrng[ Stream]  t=  Zi; 

Rand  *=  (2  *  (Zi  div  256)  +  1)  /  16777216.0; 
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end; {Band} 


procedure  AddToList  (index:  integer;  Z:  Longint;  var  Zlist:  ZstatsPtr); 
var 

tempPtr:  ZStatsPtr; 
begin {AddToList } 

New(  TempPtr) ; 

TempPtr A . Zlndex  :=  index; 

TempPtr A .  Zvalue  :=  Z; 

TempPtr  A.  Next  :=  Zlist; 

Zlist  :=  TempPtr; 
end; {AddToList} 

procedure  Randst  (Zsetx  Longint;  Stream:  INTEGER); 
begin  {Randst} 

{Set  the  current  Zrng  for  stream  Stream  to  Zset} 

Zmgf Stream]  :=  Zset; 
end;  {Randst} 


function  Randgt  (Stream:  Integer):  Longint; 
begin  {Randgt} 

{Return  the  current  Zrng  for  Stream} 

Randgt  :=  Zrng [Stream]; 
end;  {Randgt} 


function  UniformN  (x:  integer;  a,  b:  real):  real; 


var 

tempU:  real; 


begin 

tempU  :=  Rand(x); 

UniformN  :=  a  +  (b  -  a)  *  tempU; 
end;  {UniformN} 

function  Exponential  (x:  integer;  lambda:  real):  real; 

{Specify  the  number  stream  from  which  to  draw  the  seed} 

{value  and  a  real  number  for  lambda .  Returns  a  real  from  Exponential  Dist  } 
var 

tempU:  real; 

begin 

tempU  t-  Rand(x); 

Exponential  :=  -(ln(l  -  tempU)  /  lambda); 
end;  {Exponential^} 


function  gasmaX  (x:  integer):  extended; 

{This  only  accepts  an  Integer  input  but  in  order  to  allow  return  of  Numbers  for  Gammas} 
{greater  than  Gamma  of  20,  I  had  to  use  extended  numbers} 
var 
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n:  integer; 

TempProd:  extended; 

begin 

if  x  >  250  then 
begin 

writeln( 'This  is  too  large  a  value  for  a  Ganma  Calculation  '); 

end 

else 

begin 

tempProd  :=  1; 

for  n  :=  1  to  x  -  1  do 

begin 

TempProd  :=  TempProd  *  n; 

end; 

gairmaX  :=  TempProd; 

end; 

end;  {ganxnaX} 

function  GamnaN  (x:  integer;  alpha,  Betas  real):  real; 

{This  gives  a  number  from  a  Gamma  Distibution  and  is  taken  from  pg  488-89  of  L&K} 
var 

gamma,  b,  tempUl,  tempU2,  Y,  V,  Zu,  a,  q,  theta,  d,  Ws  real; 
es  extended; 
success:  boolean; 

begin 

e  :=  Exp(l); 
success  :=  false; 
if  Beta  <  0  then 
begin 

write  In  ( 'Beta  must  be  postive. ' );  {Error  check} 

end 

else 

begin 

if  alpha  <  0  then 
begin 

writeln( '  Alpha  value  must  be  postive*');  {Error  check} 

end 

else 

begin 

if  alpha  =  1  then 
begin 

gamma  :=  Exponential {x,  alpha); 

end 

else 

begin 

if  (alpha  >  0)  and  (alpha  <  1)  then 
begin 

b  :=  (e  +  alpha)  /  e; 
while  not  success  do 
begin 

tempUl  :=  Rand(x)  *  b; 

{This  P  from  pg  488  of  Law  &  Kelton} 
if  tempUl  >  1  then 
begin 
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GamrvaN 


Y  -ln((b  -  terapUl)  /  alpha); 
teinpU2  :=  Rand(x); 

if  ln(tempU2)  <=  ln(Y)  *  (alpha  -  1)  then 
begin 

gainna  :=  Y; 
success  :=  true; 
end;  {first  test} 

end 

else 

begin 

Y  *=  ln(tempUl)  /  alpha; 

Y  t=  exp(Y) ; 
tempU2  :=  Rand(x); 

if  teanpU2  <=  exp(-Y)  then 
begin 

gainna  :=  Y; 
success  :=  true; 
end;  {second  test} 

end; 

end;  {while  statement} 
end  {0<alpha<l} 

else 

begin  {alpha  >1} 

a  :=  1  /  sqrt(2  *  alpha  -  1); 

b  :=  alpha  -  ln(4); 

q  :=  alpha  +  1  /  a; 

theta  :*=  4.5; 

d  :=  1  +  ln(theta); 

while  not  success  do 

begin 

tenpUl  :=  Rand(x); 
tempU2  :=  Rand(x); 

V  :*  a  *  ln(tempUl  /  (1  -  tempUl)); 

Y  :=  alpha  *  exp(V); 

Zu  :=  sqr(tempUl)  *  tempU2; 

W  b  +  q  *  V  -  Y; 

if  W  +  d  -  theta  *  Zu  >=  0  then 

begin 

gainna  :=  Y; 
success  :=  true; 

end 

else 

begin 

if  W  >=  ln(Zu)  then 
begin 

gainna  :=  Y; 
success  :=  true; 

end; 

end; 

end ;  {while  } 
end;  {alpha  >1} 
end; 
end; 
end; 

:=  gainna  *  Beta; 
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end;  {GamnaN} 

function  Weibull  (x:  integer;  alpha,  Beta:  real):  real; 


var 

tempU:  real; 


begin 

tempU  :=  Rand(x); 

Weibull  :=  exp (In (Beta)  +  In (-In (tempU) )  /  alpha); 
end;  {Weibull} 

function  PoissonP  (x:  integer;  lambda:  real):  real; 

{Returns  a  probability  based  on  a  Poisson  distribution} 

var 

TempProd:  extended; 
n:  Integer; 

LnP:  real; 


begin 

TempProd  :=  1; 
if  x  =  0  then 

{this  prevents  some  divide  by  0  errors  -  Still  need  to  provide  some  guidance  on  input ing  } 
begin  {negative  integers} 

PoissonP  :=  exp (-lambda) ; 

end 

else 

begin 

for  n  :=  1  to  x  do 
begin 

TempProd  TempProd  *  n; 

end; 

LnP  :=  x  *  In ( lambda )  -  lambda  -  In ( TempProd ) ; 

PoissonP  :=  exp (LnP) ; 
end; 

end;  {PoissonP} 

function  PoissonNextEvent  (lambda:  real;  x:  integer):  real; 

var 

TempT:  extended; 

begin 

TempT  :=  Ln(Rand(x) )  /  (-lambda); 

PoissonNextEvent  :=  TempT; 

end; 


function  mErlang  (x,  m:  integer;  Beta:  real):  real; 
var 

tempU:  real; 
index:  integer; 


begin 

tempU  1; 
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far  index  s-  1  to  m  do 
begin 

tempU  :=  tempU  *  Rand(x); 

end; 

mErlang  :=  -(Beta  /  m)  *  ln(tempU) ; 
end;  {mErlang} 

function  NormalX  (x:  integer;  mu,  sigma:  real):  real;  {Not  a  true  Normal  but  close 
enough} 

var 

VI,  V2,  Ul,  U2,  W,  Y,  XI,  X2:  real; 

begin 

W  :=  1.5; 
while  W  >  1  do 
begin 

Ul  t*  Rand(x); 

U2  **  Rand(x); 

VI  :=  2  *  Ul  -  1; 

V2  2=  2  *  U2  -  1; 

W  i=  sqr(Vl)  +  sqr(V2 ) ; 
end; 

Y  :=  sgrt( (-2  *  ln(W))  /  W) ; 

XI  : —  VI  *  Y; 

X2  V2  *  Y; 

NormalX  2=  XI  *  sigma  +  mu;  {This  gives  Normal-' (mu,  sigma ^2)} 
end;  {NormalX} 

function  CauchyX  (x:  integer;  median,  beta 2  real):  real; 
var 

tempU:  real; 

begin 

tempU  :=  Rand(x); 

CauchyX  2*  tan  (pi  *  (tempU  -  0.5))  +  median  /  beta; 
end;  {CauchyX} 


function  fuzzyReal  (r:  REAL):  REAL;  {Only  use  this  for  a  quick  test  inAWB} 
begin 

fuzzyReal  2=  r; 
end;  {fuzzyReal} 


end . { Stochastics } 
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